#include<iostream>
using namespace std;
#include<iomanip>
#include<cmath>
#include<fstream>
#include"library.h"

double E;

void derivative(int n, double x, double *y, double *dydx){
	dydx[0] = y[1];
	dydx[1] = (x*x - 2*E) * y[0];
}

double psib(double x){

	E=x;

	int N=10000,n=2;
	double a=-5, ya=0, deriv_a=0.01, b=5;
	int i,j;
	double h = (b-a)/N;
	double *y = new double [n];
	double *ynext = new double [n];
	y[0]=ya; y[1] = deriv_a;
	for(i=0;i<N;i++){
		rk4(a+i*h, h, n, y, ynext, derivative);
		for(j=0;j<n;j++) y[j] = ynext[j];
	}
	return y[0];
}

int main(){

	double eig;
	eig = bisection(psib, 0, 1, 1E-15);
	cout<<"Ground state, numerical: E="<<setprecision(10)<<eig<<"\t theory: E=0.5\n";
	eig = bisection(psib, 1, 2, 1E-15);
	cout<<"1st excited state, numerical: E="<<eig<<"\t theory: E=1.5\n";
	eig = bisection(psib, 2, 3, 1E-15);
	cout<<"2nd excited state, numerical: E="<<eig<<"\t theory: E=2.5\n";

	return 0;
}
